

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z analyst specif


%% estimated mean utilities
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);


%% counterfactuals
load('counterfactualCH_replication.mat','srng_cfactCH')
rng(srng_cfactUS)
NS = 2e+3;
YY4 = YY(YY(:,6)==0,1:5); % 4-mpow cw-intervention strategy profiles
timeCH=0;

% CH never intervenes
VnoCH_R = UStar(YY(:,6)==0,1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VnoCH_US = UStar(YY(:,6)==0,2,1:N) + randn(size(YY4,1),NS,N);
VnoCH_UK = UStar(YY(:,6)==0,3,1:N) + randn(size(YY4,1),NS,N);
VnoCH_France = UStar(YY(:,6)==0,4,1:N) + randn(size(YY4,1),NS,N);
VnoCH_Russia = UStar(YY(:,6)==0,5,1:N) + randn(size(YY4,1),NS,N);
VnoCH_China = UStar(YY(:,6)==0,6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
noCHintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameCH',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "France" "Russia" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VnoCH_US(2:end,b,n),VnoCH_UK(2:end,b,n),VnoCH_France(2:end,b,n),VnoCH_Russia(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_France=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_Russia=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_France.^YY4(2:end,4)).*((1-sigma_France).^(1-YY4(2:end,4))).*...
                (sigma_Russia.^YY4(2:end,5)).*((1-sigma_Russia).^(1-YY4(2:end,5)));
            EV_R(m)=VnoCH_R(2:end,b,n)'*EQ(:,m);
        end
        noCHintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeCH=timeCH+toc;
end
[noCHDMES,noCHintervEq,noCHeqPayoffs] = createDMES(noCHintervEq,VnoCH_R,VnoCH_US,VnoCH_UK,VnoCH_France,VnoCH_Russia,VnoCH_China,[YY4,zeros(size(YY4,1),1)],specif);

% CH always intervenes
VallCH_R = UStar([1;find(YY(:,6)==1)],1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VallCH_US = UStar([1;find(YY(:,6)==1)],2,1:N) + randn(size(YY4,1),NS,N);
VallCH_UK = UStar([1;find(YY(:,6)==1)],3,1:N) + randn(size(YY4,1),NS,N);
VallCH_France = UStar([1;find(YY(:,6)==1)],4,1:N) + randn(size(YY4,1),NS,N);
VallCH_Russia = UStar([1;find(YY(:,6)==1)],5,1:N) + randn(size(YY4,1),NS,N);
VallCH_China = UStar([1;find(YY(:,6)==1)],6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
allCHintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameCH',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "France" "Russia" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VallCH_US(2:end,b,n),VallCH_UK(2:end,b,n),VallCH_France(2:end,b,n),VallCH_Russia(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_France=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_Russia=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_France.^YY4(2:end,4)).*((1-sigma_France).^(1-YY4(2:end,4))).*...
                (sigma_Russia.^YY4(2:end,5)).*((1-sigma_Russia).^(1-YY4(2:end,5)));
            EV_R(m)=VallCH_R(2:end,b,n)'*EQ(:,m);
        end
        allCHintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeCH=timeCH+toc;
end
[allCHDMES,allCHintervEq,allCHeqPayoffs] = createDMES(allCHintervEq,VallCH_R,VallCH_US,VallCH_UK,VallCH_France,VallCH_Russia,VallCH_China,[YY4,ones(size(YY4,1),1)],specif);


%%
clearvars -except N NS srng_cfactCH YY4 timeCH VnoCH_R VnoCH_US VnoCH_UK VnoCH_France VnoCH_Russia VnoCH_China ...
     VallCH_R VallCH_US VallCH_UK VallCH_France VallCH_Russia VallCH_China noCHintervEq noCHDMES noCHeqPayoffs allCHintervEq allCHDMES allCHeqPayoffs
save('counterfactualCH.mat')



